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ABSTRACT 


An elaust^o-pLasLic finite element, code based on small deformation 
theory is developed for plane stress, plane strain and axlsyininetric 
problems. In the plastic refflon ^^—incrementcrt flow theory Is used for 
the materials obeying Huber — von Misos yiold cr~itor-ion to model the 

stress-strain behaviour. Materials having elastlc-perfectly plastic 
and eiastlc-linearly work hardening characteristics in uni-axial 
tension are considered. Code is validated considering a V-notch stress 
concentration problem. Metal extrusion problems have been solved to 

illustrate the applicability of the code. The large deformations 
incurred in the extrusion process is modelled using Updatod Lagfongian 
formulation. Plastic >«Mrk, plastic zone spread and deformed 

configuration axM» obtained for conical die extrusion. Die surface 

friction is also considered in modelling. 
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CHAPTER 1 


INTRODUCTION 


1.1 MOTIVATION 

In i^enersl most, of t.he structural materials are Elastic- Plastic 
in nature. The complete behaviour of the materials will be well 
undex*stood if we know the behaviour of the materials in the plastic 

ranee. The non_llnaar stress-strain relations, the loailinc path 
dependence and the lar|^ deformations in the plastic ranee make the 

analysis tedious. Over the years the Finite Element Method has been 
successfully employed in anadysin^ the structural behaviour in elastic 
as well as elasto-plastic ran^^e. 

Moreover in metal forming processes like extrusion, forming and 
drawing etc., the given -material is deformed in the complete plastic 
range so as to get ttae desired shape. Visco-plasrtic models are 

commonly enq;>loyed in understanding these material pzH>cesses. In these 
processes there are two distinct regions Ca> elastic region Cb> 

elastic- plastic region. In the second region where dominantly plastic 
behaviour occurs, visco- plastic models are employed. It is necessary 
to define the boundary of this plastic region when visco-plastic 
materiad is employed in order to get the good estimate of power 
required, finished pxoduct shape and forming pressure. 

Keeping the above points in view the present thesis is mainly 
devoted to the elasto-plastic finite element implementation and 
extending towards the application of nmtal extrusion problems. 
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±J2 Aim And Scope Of The Peeacent. Work 

An increments FTnlt^e Element.- An^ysis has been developed for t.he 
elast.o-plast.lG problems. The plasblc st.ress-st.raln relations In t.he 
plasbic ranipe are modelled uslne Pr'andtl~R»uss relation. The following 
methods are en^loyed In solving t.he final dlscx^t.ised conblnuum 
equations viz. NGUfton-Raphson method, HcKHfi^d N»wton~Raphson method 
and Yamada's EU approach. In t.he later ai^roach an increment load 
sufficient enou{;h to yield one element at a time is applied and the 
loads ai*e inczemented accordingly. Ck:>nstant stress/straln triangular 
<3 - nodal> elements have been employed. Materials having 

characteristics in uni-axlal tension as elastlc-perfectly plastic and 
elastlc-llnearly work hardening are considered. 

1^ Outllxke Of The Thesis 

In the present chapter motivation, aim and scope and out line of 
the thesis are presented. The review of the literature for the present 
work is given in cha 4 >ter 2. In chapter 3 the elastic and plastic 
stress-strain matrix, finite element formulation for olasto-plastic 
analysis and all consrtitutive relations for plane-strain, plane-stress 
and axlsymmetric cases aure derived. Some of the numerical methods for 
solving non-linear problems are presented in the same chapter. 
Formulation and flow sequence of the program for the extrusion 
problems axe given in chapter 4. Results and discussion axe presented 
in chapter 5. In chapter 6 conclusions and scope for the future work 
axe presented. 
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CHAPTER 2 


LITERATURE SURVEY 


Har-cal and King 141 su^^est^ a slirfness met^hod for 
elasi.ic-plast,ic problems. They obtained the Incremental sttess strain 
relations for an elastic-plastic materials in terms of par-tial 
stiffnass coefficients. The partial stiffness coefficients for the 
solid body of revolution with symmetric loading are ^Iven as follows. 

With cylindrical coordinates r> z and usii^ the conventional 

notations _ for the stresses, the von Mises yield criterion can be 
written as 141 

Co - o >* + Co - o >* + Co - o >* ■ 2 o * C2.1) 

Or r z X O 

Equation C2.1> can be put in an implicit differential form as 

3 o' do + 3 o' do + 3 o' do + 6 o do ■ 2 o do C2.2> 

O & rr zz rzrz 

where 

o' denotes the deviatoric stress e.c. o' ■ — C 2o -o -o > 

0 a 0 r 2 

o is the equivalent stress 

The prefix d in eq. C2.2> is used to denote the changes over the 
small Increment of the load. 

Marcal and Kin^ C4] used the foUowini; incremental stress strain 
relation 

C<h:> - to 3 <do> C2.3> 

m 

where Is the 5X5 symmetric matrix. <do> and Cdc> are 

column vectors of stress and strain. 
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Equat^ion C2.3> can be written as 

<do'> ■ to- l~*<dc> C2.4> 

TO 

[O' 1 ^ is called as par-tial stiffness matrix and cam be wrlbben as 

Tf\ 

da, 

m — ~ Tor i>J ■ e, r, z, rz^ p 

j 

Here a is used t,o represent, s ^ and c equal to zero. The 
P p 

coefficlen'ts / dc ^ ai'e called partial stiffnass coejffitriants. 

In usin^ bhe above paorbiad. sbiff ness mebhod^ complications arise 
when the elastic part of the structure with stresses near the yield 
stress becomes plastic with the next increment of the loaul. The re^^ion 
which is adjausent to the elastic- plastic inter fau^e had been called the 
transition ra^on. For this Marcad and Kin^ C43 defined a mean 
stiffness of the element by wei^htin^ of the elastic and 
elastic- plastic stiffness in the foUowini; ratio. 

Co ■ m Co 3* + <l-m> Co 3*^^ 

TO TO TO 

where the superscript ep indicates elastic-plastic and m is the 
strain required to yield / estimated strain for the inclement of 
loaKl 

Later Yamada, Yashinutra and Suktjtrcri C13 obtained an explicit 
relation for CD**^3 <i.e. elasto-plastic stress strain matrix> for von 

Mises materials usin^ Prandtl-Rauss relations. The expression obtained 
takes a simple form and can be readily accommodated to the finite 
element anadysis. At each Incremental staq^e of cadculation which 
traces the expansion of the elastic- plastic interfatce. This method is 
similar to that of Marcal and Kingi41, but uses the smadl and varying 
increments of load is applied in every iteration. 

The procedures of cadculation and flow sequence of Yaunaula's 
method are i^iven in chapter 3. The present thesis uses the Yamada's 



method for elas to- plastic analysis. 

Lea and Kahayashx IZ\ analyse the flat punch indentation problem 
for elastic- plastic boundaries and for stress distributions by usin^ 
Yamada^s algorithm. 

Osakad and Fujina C33 did the Elas to- Plastic analysis of 
hydrostatic extrusion usin^ the same Yamada^s technique. 



CHAPTER 3 


FORMULATION OF ELASTO-PLASTIC PROBLEM 


In -this ch3pt.er necessary const-1-tutlve rela'tlons for ma-terials 
obeying; Huh^r^van ^#is*es yielcf cr^it^r-xc^n and Tollowln^ Firandtl-R^uss 
xncr'&m&ntal /lou> theoryCJ^-incremenbal plasbldby flow t^heory/ Von 
Mises riow rude > are presented The maberiai is assumed to be 
homogeneous and isobropic during the entire loading path. The behavior 
of the material in uniaxial tension is characterized by Elastic - 
perrectly plastic and Elastic -linear work- hardening nature. 

Fi^. 3.1a shows the uniaxial tension stress-strain relations Tor 
Elastic-perfectly plastic model and are expressed mathematically as 


s ■ 


a 

E 


for a < a 

y 


C3.1> 


C m 


E 


X 


for 0^0 

y 


C3.2> 


where 


E 

s 

a 

a 

y 

X 

Fli?. 3.1b 
work hardening 


is younc^s modulus 
is uniaxial strain 
is uniaxial stress 
is yield stress 

is a scalar to be determined and is greater 
shows the stress strain relations f or 
model and are expressed as 


than 0. 
elastic- linear 



o 

y 

E 



-for a < a 

y 


for a > a 

y 


<3.3> 


■» 


C3.4> 






whe]>e 


is tar^enti^ modulus. 

The Small strain th»ar^y Is applied for every increments load. 
The loading are monobonic and increasing in nature. The Loadings are 
also considered to be proportional in nature. Updatad La^ang±an 
farmulatian is used to ^et the deformed configuration. 

In the following sections . represents tensoriad strain and <£> 
represents en^^ineerin^ strain vector. 


s «-<u + u. .> 

1-J 2 v,J JA 




xy 


yz 


V J'x* 


du 

dx 

an 

y 

ay 

au 

2 


au 

> 

dx 

au 

> 

az 

au 

2 

a5r 


au 

X 

a^ 

au 

2 

dy" 

du 

> 


where u , u and u are ftinctlons in x> y and 2 which |dve 
X y z 


displacements in x^ y and z diiections respectively. 


3.1 Elastic Stress-Strain Matrix 

The total strain increments are made up aT increments 

d£:*.of the elastic strain components auid incieinents ds^. of the 

VJ 
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plastic strain components. 

d£ . - dc* .+ d£^ C3.5> 

i- j vj Lj 

The elastic strain Increments are related to the Increments of" 
stress <da > by Hook's law. 

■ = — do*. . -* — do, , 6 C3.6> 

Jb vj b kk vj 

The eqxjatlonC3.6> can be solved In terms of stress increments 


at 

Lj a+v> 


+ 

'-j 


vE 




<3.7) 


Strain increments can also be written in terms of deviatoric 


Co-' .>stress. Equation C3.6> becomes 


J 



1-t-u 

E 


do' 
v j 


Cl-2u> 

3E 


kk \.j 


C3.8> 


Where 

u is the poisson's ratio 

f 1 if i ■ j 

6. . is Kronecker's Delta « -I 

t 0 if i ffc J 

The equation C3.7^ can be written in the matrix form as 

<o> ■ CD®1 <£> <3.9> 

Where 


CD^] is the elastic constitutive or elastic moduli matrix 

<£> is the elastic strain vector 

T 

<£> m £ £: s £ ^ y y 1 

X y z xy yz zx 

<o*> is the stress vector 

T 

<o> »Cooo' T T rl 

X y z xy yz zx 
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CD"3 


Cl+v>a-2v> 


Cl-v> 




Cl-v> 


Symmetric 


O 


Cl-2v> 


0 


0 0 

Cl-2v3 


O 

Cl-2v> 


Equation C3.6> for total stress and strain can be written in the 
matrix form as 


<£*> ■ CG*3 <o> 


O.IO 


Where 


CG*] is elastic 


compllence 


mat^rlx 


cc*] . 4- 


1 -V 

1 

Symmetric 


-V 0 

-V 0 

1 0 

2Cl+i^> 


0 0 

0 0 

0 O 

0 0 

2<l+i>> O 


2<l+y> 
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3«2 Eiasto*Piast»iG StreaRs-SIrain matrix 


An axpUcit axprassion of alas tic* plastic strass*strain 
matrix CD^^l had been obtained by yomadcTCll The plastic strain 
increment in terms of plastic potential function ^Co. .> is 

l^iven by C13 

ds^. - ^ 

I-J ^cr. . 

dX. is a positive scalar f actor of proportionality^ which is non zero 
only when plastic deformation occurs or plastic loading occurs. 

The associated flow rule is obtained by considering yield 

funGtlon<f> for plastic potential in equation C3.11> 


■ dX — 

V j aa. 


C3.12) 




fto" . > - J - 

VJ ** VJ ''z 


whexe 


K 


is the yielding shear stress 


J is second invariant of the deviatoric stress tensor 
Equation <3.12> takes the simple form as 






dX O'' 


<3.13> 


The above equation is called Fr^andtl-R^uss equation. 

Using equations C3.5> and C3.8> the Prandtl-reuss equation for 
the total deviatoric strain increment d^' . during continued loading 




are expressed as 


/ do. 

dc' mcr. .dX+ 

"J 2G 


<3.14> 


Where 


dX 


3 

2 


de 


P 


a 

The derlva*tion for dX is ^Iven In 
Equivalents st^ress a and piastsic strain 


3 da 

2 - H' 

a 

tshe appendix 
increments ds 


-I. 

^ are given by 


<3.15> 




a' a' 

^ 3 


C3.16a> 


- 2 
a 


i[ 


Ca - a )^+ <a - a !>^+ Ca - a 2 t ^ + 2 t ^ + 2 t ^ 

X m y m z m xy yz zx 


Where 


a + a t a 
X y z 


Afber Simplif icabion 


2 2 2 

a + a + a - a a - a a 
X y z X y y z 


- O’ O' + 31 ' 
z X I, 


2 ^ 2 ^ 2 
T + T + T 
xy yz zx 


H - 


da 


ds 


/ 3 IJ LJ 


<:3.16b> 


corresponds bo bhe slope of the curve equivalents 


d£ 


sbress a Vs plasblc sbrain CSds do* is bhe increments In tshe 


equivalents stsress 

The Von Mises yield crltseria is given by 


a, , a . 

V J V J 


2-2 
3 ^ 


Taking; bhe dlfferen-Lials on both 
equation 


<3.17> 


the sides of the above 
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cr' do-' • § O' do- 

i j I j 3 


Making use or equa-tion C3.15> rhe above equation can 


O'' do-' - J ^ * H dk 

'■ i >■ J V 


Eliminating do' . from equa-tion <[3.14> and equa-tion C3.18) 


'• J 


2G o- 


' (de' . - o' d\ 1 ■ ~ 

«-jL i-J >-j J 9 


o * H dX 


Rearranging -the -terms to solve for dX 


dX 


[I 


o * H + 2Q 


cr' . a' . \ 
i-j v J J 


28 o' de' 
>■0 >-J 


2G o' de' 

i- j 


dX 


^ o * h’ + 28 o' . o' . 

9 vj vj _ 


dX 


O'' . dc' . 

VJ V J 


2-2 

3 


i -4- 

H' 


36 


dX 


o' de' 

i-J 


With 


dX 


^ 2 - 2 r 

^ Mi —.1. f'y 1 

1 + - 

H' 1 

^ 3 ^ [ 


38 J 

can be wrib-ten 

as 


O'' de' 

^ j _ 

<y* 

T- j 

de . 


o' de' « o' de since o' . » o' + o' + o' is identically 


be written 


C3.18> 


<3.19) 


<3.20) 


<3.21) 

zero. 
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<:3.22> 


Devlatorlc strain tensor d^:' . is defined as 

^ j 

d^: 


d£ ' ■ ds . - S - ^ 

1- J l-J d 


kk 


Where 

Hx: m + d£ + ds 

kk X y z 

Now substituting dX of eqpjation C3^1> back into equation <3.14> 
and making use of ecjuation C3-22>^ the deviatorlc stress increment 
tensor can be expressed as 


do* « 2G 

j 


1 fir' - - 

Im Im 

1 Cf . , 

t 1 j ^ J 

s J 


ds. 


do* . ■ 2G Ids: - <5 - 

^ J I ''J 3 


kk 


o* . o* de 

V j L Tn Itn 


<3^> 


The -to-tal st,ress incremen't da. . is defined as 

'■j 


da. 


'-i 


da‘ + 

^3 


-ii da 
3 kk 


do-' + 
^3 


3a-2v> ^ii. 


^ ^ 2a+v>G , . 

da‘ . + — — 6. de. 

vj 3Cl~2v) vj vv 


Finally, subs-ti-tu-ting for- do-' . from equation C3^3> yields 


J 


do ■ 2Q dc. . + 6. . de, . 

1 . j L ‘•J <l-2v> VJ VJ 


‘^kl “^kl 1 

" S J 


C3^4> 


The equation is rapresen-ted in mafrix form as 


<do> ■ CD ^3 <de> 


C3^5> 
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The plas'tlc st.ress-s'traln ma'trlx CD ^^3 is symme-tric and is 
expressed in t.he explicit, rorm. 
tD^^l is ^iven by 



Where 

a-i>> 

^ “ <:i-2v>> 

H ■ ^ 

a-2v> 

3.3 Finit.e Element. Formulet.ion 

The general governing equat.lon of the finite element method for 
static analysis can be derived from the principle of xji.r-tttal wor-h. 


I O' Se dV ■ I F 6u dA + I q.6u,dV 

J ij Ij J " J " 

V A V 


C3.26> 





where 6u 


vlrt^ixal 


strain 


and 6c . are virtual displacement and 

L J 

increments, and they form a compatible set of deformations. F. 

i, 

and q are surface traction and body forces,V and A represents volume 

and area of Integration respectively and a with F and q form 

i j L L 

the equilibrium set. In a matrix form Eq.C3.26> becomes 

<6u>’^<F> dA + j <6u>’^<q> dV C3J27> 

V A V 

where the vectors for displacements < u > , strain < ^ >, and stress 

<€/> are defined as 


T 

<V1> 


[ 

U 

u 

u ], 

<6u>^* C 6u <5u 6u ] 

C3^> 




1 

2 

3 

12 3 


T 

<£> 

■1 

[ 

£ 

C 

c 

y y y ] 

C3^9a) 




X 

y 

z 

X y y z z X 


T 

<6£> 


c 

6c 

6c 

6c 

6y 6y Sy 3 

C3^b) 




X 

y 

z 

xy yz zx 


T 

•Ccr> 


[ 

a 

a 

a 

T T T 3 

<3.30> 




X 

y 

z 

xy yz zx 



For a geometrically linear auialysis or a small def ormation 
analysis, we have, 

<s> ■ CB] <U> <6£> ■ CBI <<5U> <3.31> 

where <U> is the displacement vector for nodal points that is related 
to the distributed displacement <u> by, 

<u> - CNl <U> <3.32:> 

in which C N 3 is the matrix of clx:splac&m&nt int^irpolatxon function,. 
or the shape /unction > and the strain displacement matri,y: CBI Is a 
matrix defined as 

CB] » CL] CN] c:3.33> 
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iao^d CL3 Is the dx//erentiaZ operator matrix dertned ats 


L - 


dx 

0 

0 

dy 

0 

dz 


0 

ay 

0 

a 

ax 

a 

az 

0 


0 

0 

a 

az 

0 

a 

ay 

a 

ax 


such that 


<£> * CLl <u> 

5^abstitutin^ the Eqs. <3,31> and C3.32> in Eq. C3^7>^ 
equation for a small- deformation analysis is 

I dV ■ [^n«’*'<F> dA + 

“A 

or 

CB]’'<<y> dV ■ <R> 

■^v 


CN]’^<q> dV 

V 


where <R> is equivalent, ext,ernal foi>ce acting on the nodal 


<R> 


r CN]'*^<F> dA + r Cro’*‘<q> dV 
"^A V 


<3.34 > 


<3.35 > 


the governing 


<3.36a> 


<3.36b> 


pyolnts. 


<3.37> 
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Therefore ^usln^ < 0 '> ■ CD]<£r> ■ CD3CB1<U> "the governing equ 3 t>ion ceoi 

be written as:> 

CKJ <U> * <R> C3.38a> 

In Incremental form 

CkJ <dLr> - <clR> <3.38b> 

where CKD Is the stiffness matrix of the structure , 

era - I CBJ'^CD] CB] dV <3.38c> 

in which CD3 is the constitutive matrix. 

For elastic case CDl is replaced by CD^I and for plastic case CDl 
is replaced by CD^^l. 

The whole domain is divided into triangular elements. The 
following figure shows one such triangular element with three 

nodes i, J and m. 


J 

The inat.rix CB] Is given by 

CB] ■ CB B B ] 

i j m 
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0 


- -2X 


Where 




X - X 

m j 


X > 

V 

elemen't 


th 


y are "the x and y co-ordina'tes of "the i iKKie in t.he 


A * Area of the the element 


-fcxy -xy> 
2 2*^1 


< X y 


X y > +C X y 

9 2 9^1 


X y > 

19 


Incremental Analysis 

I 

In this analysis equiation <3.38b> has been converted into 
era <AU> - <AR> 

cia"’< - u™ > ■ < - r'"> 

era is evaluated at increment and the problem is solved for 

Um+i 

In an incremental analysis^ the total load <R> acting on a 
structure is added in Increments step by step. At the <m+l>^^ step^ 
the load can be expressed as 


Tn+ 1 


<R> 


"<R> 


Tn+ 1 


<AR> 


C3.39) 


where t,he left, superscript. m has been used t.o Indicate bhe 
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incremental step. Assximln^ the solution at the step^ 

are known^ and at the Cm+1>^^ step^ corresponding to the 
load increment <AR>> 


U > - U > + < AU > 


C3.40> 


a y m a > + < Aa y 


C3.41> 


Here, the left superscript for the increments has been dropped. 
Eq. <3.36> becomes 


m+ 1 


< F > 


m+ 1 


< R > 


C3.42a> 


where p > equivalent force of stress acting on the nodal 

points, 

-f CBl’^ <0'> dV <3.42b5 

•^v 

or 


I 


[B3 <Ao-> dV 


m+1 


<R> 


-I 


[sf "'<C7> dV 


<3.42c> 

V V 

Eq. C3.42c>, In fact, represents the equilibrium of the external 

forces, with the internal forces Two types of 

algorithms are therefore involved in solving Eq. C3.35> for the stress 

increment < A<7 > and displacement increment < AU >. 


3.4 Special Cases 
1. Plane Stress 

Elastic case 

In this case 


20 



<y = 0 ; T = 0;t = 0 and 

z yz xz 

r • 0 } r “ 0 ; ^ C cr + a > 

yz xz z L X y 

The stress vector Is reduced to 

<o> * C O O’ T 3^ 

X y X y 

The strain vector is induced to 

<£> ■ C £ £ 3^ 

X y xy 

Elastic moduli matrix 



Elastic compliance matrix 
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Whkure 


2 H' “ 2 ^ 

O' + 


xy 


PE 


Cl+v) 


R ■ a' ^ + 2vo'' cy' + cr' ^ 

X X y y 


Q ■ R + 2P<l-v > 


O '* (< 7 ^ + 0^-00 + 3 t ^ 12 

lx y X y xy J 


i{ 


2 2 2 
S ■i’ S •i' £ - £ £ 

X y z X y 


z y X z 


4 } 


2.Pl3ne S-train 


Elaartic case 


In "this case 

£ m 0 Y 
z yz 


0 ; r 


T ■0;t “0; am 

yz zx z 

The st^ress vecter Is reduced *to 


0 ; and 

E i> 


a+v)<i- 2 v> 


£ > 
y 


<0> * [ O O T ] 

‘ X y xy 

The s-train vector is reduced to 


<£>■[£ £ r 3 

X y xy 
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rHiCSJ 



ED*] 


<l+v>Cl-2v> 


Cl-v> 


a-v> 




EC *3 


a+v> 


Cl-v> 


- 1 ^ 


<.i-uy 

0 


For pl£tst/lc C£ise 


ED*»^3 


0 

O 

2 





o' o' 

X y 

o' T 

_ X xy 


<:i-2v) s 

<l-2v> 

s 

S 

E 


Cl-v> 

o'^ 

y 

o' T 

_ y xy 

Cl+v> 


a-2v> 

S 

S 


Symmetric 



2 

. T 

1 X y 


In this case stress component in z-direction is not equal to zero 
and is given by 
E 


do 


[(<!-: 


a+v> l<i-2v> 


- ^ M 


a a' ^ <y' T 

zr> - 


Eqmvalen-t s-tress a is ^iven by 


2 2 

O' + O' 

X y 


o - o o - a a 

z X y y z 


^2 J 2 

O O + 3 T I 

z X xy J 
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and £ is given by 


2 r 2 

3 [ 


+ £ - £ £ + 
y X y 


3 1 - 

i ^xy J ^ 


[B3 ma-trix for fhls case Is same as in -the plane sfress case 

3 . Axlsy iime-tr y 

Ellas'tic case 
In -this case 

’ ■"rS ■ “ > 

■ ® = ""re “ ® > 

The sfress vecfor becomes 


<ay ■ C <y o O' T 3 

r z e rz 


and bhe sfraln vector becomes 




< e c y y 

r z e rz 


r and e represent . the radial and tangential 
respectively. 


directions 
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[0*3 


a+v><i-2v) 




V 


Cl-v> 


Symme'tric 


V 


<l-v> 


0 

Cl-2v> 


EG 3 ■ 


1 

E 




Symmet.rtc 


For plast^ic case 


-t> 


-V 


O 

O 

0 

2a+v> 


[0*^3 


Cl+i^> 


<l-v > 
<l-2v> 


Symme'tric 


cr a 
r z 


Cl-2v>> £ 

* 

Cl-v> ! 

Cl-2v>“ S 


r e 


Cl-2v> 

V 

ci-2i^y 

Cl-v> 


2 e 


s 

2 


O'' T 
r rz 


O" T 
2 rz 


e rz 


1 

2 


T T 
rz rz 


3.5 NUMERICAL ALGORITHMS FOR SOLVING NONLINEAR EQUATIONS 


In -this section ,-the two methods of the Newton type that have 
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been widely used In finibe element analysis are discussed. 


Considering that the stress < o* > is a nonlinear function of 

displacement y <IJ>^ Eq C3.42a> can be rewritten as 

U > ) - F < U > > > - R > C3.44> 

Eq <3 .44!) is a nonlinear matrix equation , expressed in terms of 
displacements y y ^ noted previously in Section 3.3^ this 

equation represents an equilibrium of the external forces, 
with the internal forces , The iterative method used to solve 

Eq. C3.402> is therefore called the eqfuilihrium iterotit^e method. 


3.5.1 


Newton-RapHson Method 


Ve have already obtadned the approximation, ^ ^<U>^ ^ 

to the displacement Expanding usin^ the Taylor 

series expansion at ^ ^^and ne^lectin^ all hig-her order terms, 

we get 




± ) 


> + 


a* 

du 




ir>+l_aj>< v-l ) ^ 


■ TT1+ 1 < i - 1 ) 

<U> 


or 




<u> 


< i > 




m+ 1 


<R> 


' m+ 1 < t - 1 ) 

<U> 


where 


0 


C3.45> 


<3.46 > 
<3.47> 
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Recognizing that 


« t “ 1 > 


dF I 

du 

' m+ 1 < i. - 1 ) 

<U> 

- I CBl'^ CD*^]| [B] dV 


C3.48> 


where < i. - 1 > the elastic- plastic stiffness corresponding 

^ <u> 

1,0 displacement. and tan^ent,ial 

st.irfness matrix of the structur*e^ ve obtain the iteration scheme of 
the Newton- Raphson algorithm as 





m 


+ -CAU>‘ ' 


m 

”'-cu>. 



m 

"'<F> 

a ■ 


C3.49a> 

<3.49b> 

C3.49c> 


This Iteration continues until the proper convergence criterion 
is satisfied. The iteration procedure is shown schematically in 
Fig. 3.2, 

The Newton- Raphson method converges quadratically at a faster 

rate. However, it should be noted from Eq. C3.49> that the tangential 

stiffness matrix , is evaluated and factorized at each 

iteration step , and such an operation can be prohibitively expensive 
when a large-scale system is considered. Moreover, for perfectly 

plastic or a strain softening material, the tangential stiffness 

matrix may become singular or ill-conditioned. This may cause 
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difficulty in the iteration procedure. Modifications of the 

Newton- Raphson algorithm are therefore necessary, and will be 

described in the following section. 

3.5.2 Modified Newton- Rap lison Method 


One of the modification of the Newton-Raphson method is to 
replace the tangential* stiffness matrix K in Eq. <3.49a> 

by K ], which is a tangential stiffness matrix evaluated at load 

step n < m+1. Matrix K 3 is evaluated at the beginning of each 

load step, or for the step, the stiffness matrix 

- '\ki <;3.5o> 

is used. The iteration scheme for the modi/ied N^tM^an-RapHsan 
algorithm is expressed as 



<3.51a> 

C3.51b> 

C3.51c> 


Again, the Iteration continues until a proper convergence 
criterion is satisfied. This modified iteration pin>cedure is shown 
schematically in Fig. 3.3. 

The modified Newton-Raphson algorithm involves fewer stiffness 
matrix evaluation and factorization steps. As a result, computational 
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Fig. 3.2 Newton-Raphson Method 


F 



Fig. 3.3 Modified Newton-Raphson Method ' 



effort for one Iteration cycle is much less than for the 

Newton-Raphson algorithm for a lar^e-scale system. However, the 

modified Newton- Raphson method converges linearly and, in more 

slowly than the Newton-Raphson method, l.e., for a ^iven problem, more 
iterations are needed to reach a convergence when the modified 
Newton-Raphson method is used. The convergence rate of the modified 

Newton -Raphson algorithm depends to a large extent on the number of 
times the stiffness matrix is updated. 

Moreover, the problem that the stiffness matrix may be singular 

or ill-conditioned still exists. Another problem associated with the 
modified Newton-Raphson method is that if a change in the external 

load cause an unloading, i.e., the stress state is unloaded from 

plastic state to elastic state, this method may not lead to converge 

for that iteration , unless the structural stiffness matrix is updated 

for this situation. This in turn increases the complexity of coding in 

the numerical Implementation of the modified method. 

Convergence Criterion 

A properly defined convergence criterion to terminate the 
equilibrium iteration is an essential part of an efficient incremental 
solution strategy. At the end of each iteration , the solution 

obtained must be checked against a selected tolerance to see whether 
convergence has occurred. 

For aui equilibrium iteration solution < U > is sought at which 
the difference between the external force and Internal force, 
R > - F > >^‘*less than or equal to an allowable value. This 

criterion is called /or-ce criterion. 
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The -tolerance must be carefully chosen. A loose tolerance ^Ives 
inaccurate results, while a too tight tolerarice may leads to wasteful 
computation to obtain a needless accuracy. 

3.5.3 YAMADA'S METHOD 

In this section a step by step incremental algorithm to solve 
non~linear problems is presented. The sec[uence of the program is as 
follows 


1> Apply a unit test load or unit displacements and compute 
the displacement vectorCAu>, <y* and <y 


2> Find 

bhe 

elemenb 

which 

has maximum oCsay for elemenb 

>. Let 

a be 

bhe 

maximum 

a of 

all bhe elemenbs. calculabe a scale 

factor 

max 






according 

bo 

r ■ Y/ 
© 

a 

max 




3> Multiply Au, e, a, a and load <AF> by the scale factor r to 
ob-tain bhe values corresponding bo bhe load bhat is required bo Jusb 
yield one elemenb. The elemenb <i.e. bhe elemenb which is having^ 

maximum sbi^ss> is considered as a posb yielded elemenb from bhe nexb 
compubabion. 

4> Add Au bo bhe inibial coordinabes of bhe nodal poinbs in order bo 
follow bhe pabh of deformabion. 

5> Galculabe elasbic^platsblc sbiffness mabrix for all bhe 

yielded elemenbs using CD* ^3 given in secbion 3^. Replace bhe elasbic 
sbiffness mabrlces wlbh bhe above calculabed CK 3 for bhe yielded 
elemenbs in bhe global sbiffness mabrix. 

6> Apply an Incremenbal load or displacemenb and compube bhe 
corresponding values of Au> Lb and La. 

7> By making use of bhe following equabion calculabe r for every 
elemenb remaining in bhe elasbic sbabe. 





C3.52> 


r + -/ r 4 <Aa 

8 t op 

r * 

2CAa 

O t 8p 

Derivat^lon Tor r is given in the appendix -II 
8> Find the element which is having minimum r Csay Tor element > 

and let be the r f on -that element. The load incnement multiplied by 
the factor is just sufficient to yield the element. From the 

next iteration element is considered to be yielded. 

9> Multiply Au, Ao* and load increment by the factor r and a dd 

to the present displacements > strains > stresses and load. Add Au to 
the coordinates of the nodal points. 

10> Compute a for each element. 

11> Check a of any element remaining in the elastic state has reached 
0.995Y. The element having a is equal to or greater than 0.995Y is 
treated as yielded element in the next step and thereafter. In doing 
so total number of Iterations will be reduced. 

12> Calculate the increment Ai:^of equivalent plastic strain for 
each post yielded element using 


a2 ^ 


a' As 

^ j 



r H' 1 

O' 

1 

— — 1 

1 CO 

+ 

T-l 


13> Check that As^ is positive for all the post yielded elements. If 
positive^ go to the step <5> and rep>eat the process otherwise stop 
doing the process. 


As^< 0 means the plastic region of the non hardening material 
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expands t^o such an extent/ t/hat. large plastic strains are possible 
without any practical increase of load. Though the implementation of 
the ^ has not been f ully understood, it can be concluded that it is 

aoi indication that we approach very close to the collapse load. 
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CHAPTER 4 


FORMULATION OF THE EXTRUSION PROBLEM 

This chapiter deals vith Elastic-Plastic analysis of extrvislon 
process .The following -two cases aio dealt. wit.h 

1. Transient, analysis. 

2. St.at.ic analysis. 

In transient, analysis the die is not completely filled with 
material as opposed to the static analysis in which the die is 
completely filled with the material before starting the process. These 
cases are analysed for both plane-strain and axLsymmetrlc problems. 

The flow sequence of the program and the expressions for total 
work done, plastic work done, frictional losses and the normal die 
reaction forces are discussed in detalL 

The total work done is calculated by using 

s 

t ^ J 

W ■ / / o- dc. . dV 

V o 

Where 

V is 'the volume or the domain 

d£ m d£* + de^ <4.2> 

ij i-j '■J 

The els^rtlc work done is given by 

. 

W* ■ f X O' de. , dV C4.3> 

V o 

The plastic work done is the difference of the total work 


the following relation 

<4.1> 
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done and elast^ic work done 

Piast^lc work done and elas'tic work done are shown in figure 4.1 
at a particular stafe of st-ress. 

41- TRANSIENT STUDY 

The figure 4.2 shows fhe billet and the punch positions. The 
billet is forged at the end to fit into the die before the extrusion 
process takes place. The discretization of the material is shown in 
fi^, 4.3. The symmetry about x-axis is exploited to reduce the 
computational cost. A finer mesh is taken near the die surface to get 
good results.The domain is discretized into triangular elements. 

A computational algorithm for the incremental displacement method 
using the elasto- plastic analysis for the above problem is presented 
here. 

1> Apply a unit punch displacement and compute elastic 
displacements at the nodes ‘Cu> , forces at the nodes strains 

and stresses <o’> at each element. 

2> Calculate the equivalent stresses a and equivalent strains ^ 
of each element. 

3> Find the element which has the maximum o* Csay Tor J 

elements. Let <y be the maximum o* of all the elements. Calculate a 
max 

scale factor r ■ — — ^where Y denotes the yield stress of the 

O “ 

a 

max 

materiad under consideration. 

The folio wine figure shows the factor r^ graphically. 
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Fig. 4.1 El,astic Work and Plastic Work 
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4. S The pjosition of material, punch and die when 
the die is not completely fill Ctransient case3 






str&ih 

4> Multiply <u>, <F>, <<y>, <£>, a, and c by the scale factor r 

to obtain the valvies corresponding to those at the elastic limits. The 

» 

elemen-t <'the element which is having a > is treated as an 

TTiax 

yielded element for the next computation. 

5> Add <u> to the initial co-ordinates of the nodal points 
in order to follow the path of deformation. 

6> Calculate CD*^] and hence forth CK*^] for all the 

post-yielded elements or continued loaded elements [D*^3 is given in 
chapter 3. Replace the stiffness matrix for the yielded elements with 
the above corresponding plastic stiffness matrix. 

7> Apply the Incremental punch displacement and compute the 
corresponding values of fAu>^ <AF>> <A£>> <Ao'>. 

8> By making use of the following equation calculate r for all 
the elements those are remaining in the elastic state. 


r + y r 4 <Ao <y^ - 

st«p 


2CAo > 

a t ap 


Find the minimum r Csay for element r is mlnimum>. Let the 

minimum r be r 

1 
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9> Calculate Tor the node which Is nearer <say 1^^ node> to 

the entrance of the die and on the container surface using the 
following relation. 

< X - X > 

dtdont V 

r « 7 

2 Au 

XT. 

Where 

x_, ^ is t.he x-coordinate of the die entrance. 

is the X” coordinate of the node which is nearer to the 

die entrance. 

is the displacement in x-direction of i’^'^ due to the 
Incremental punch displacement. 

r^ determines the punch displacement that is Just sufficient 
for the i*'*^ node to reach the die entrance. 

The following figure shows the notations used to calculate r^ 





10> Calculate r of all the possible nodes which may come in 

C 

contact with the die surface due to the incremental punch displacement 
by making used{the following relatlorvs. 
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Y 

1 


a 


b 


* Y 


d\.« e n 1 


Y « Y 

2 dleent 


Y - 
1 




Y + Au 
j yj 


t.an<a> C x - 

j 

t.anCo() C 


Y 


2 


X > 

dtttonl 


X + du 

J xj 


C4-50 

.. > , (4-5b) 

d T. ©©n t 


r 

CJ 


a 

a + b 


C4-^e; 


r indicates "the factor r 

cj C 

Where 


of the 


node. 


Y_, _ is the inlet radius for axisymmetric case or half of 

d t « e n t 

the inlet thickness of the billet for plane-strain problem. 

The following figure shows the all notations used to calculate 


r . 

CJ 




Let the mlnlmumCr ) be r <say r is minimum for k node!). 

cj 3 C ■ 

11> Let the mlnlmxjun of r , r and r be r . 

12 3 mvn 

12> Multiply <Axi>, <AF>^ <Ac?'> and <A£^> with the factor 
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13> Add 'CAu)', ■CAF!)', ■CAo’!>, ■CAc)- t^o t.he previous st/at.e {uJ, 

<C'>, <£> 

14> Updat-e t-he co-ordinates of all the nodes. 

15> Calculate the equivalent stress and the equivalent strains of 
all the elements. 

Calculate the Incremental total work done dW*^ ^ elastic work 
done V . Add dV to V to get the total W , and then calculate 
using the relations given at the beginning of this chapter. 

17> conditions to be imposed from the next iteration 
<1). If r is r then^ 

1 mtn 

element^Ci.e. t^he element for which r is minimum^ is 

considered 3s the post- yielded element from the next computation 
<2>. If r is r then^ 

2 mtn 

node Cwhich is nearer to the die entrance end on 

the container surfaGe> is reached the die entrance with t-his punch 

displacement. From the next Iteration node is considered to move 

alone surface CD. 

O). If r is r . then^ 

3 mtn 

node Cwhich is havine minimum r !> will come in 
contact with this pxinch displacement. From the next Iteration k 
node is considered to move adone the die surface CD. 

18> Check if t-he displacement, of t-he pvmch is reached bhe 

specified displacement, t-hen shop doing t-he process^ ot-herwise go bo 

bhe sbep <6> and conbinue bhe process bill bhe above condlbion is 

sabisfied. 

BOUNDARY CONDITIONS; 


The boundary condibions for bobh bhe plane sbraln and axisymmebric 
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cases 


alon^ t/he surfaces AB> BG^ CD> DE^ EA are assumed fo be as 


follows in regard bo bhe nodal forces F^> and the nodal 

displacements u^>. 

AB; u * unit punch displacement and F *0; 

X y 

BC: F ** 0 and u * 0 ; 

X y 

CD: displacement normal to the die surface is zero 


l.e. u tanCot> + u ■ 0 

X y 

and 


We know that F * /l/F 

l 

F * F cosCa> - F sinCot> 

Lx y 

F ■ F sinCo(> + F cosCa> 
n X y 

F C utanCot> - 1> + F <. ^ + t.anCa> ) ■ 0 

X y 

It is assumed that along the die surface there is no 


sticking and only slip occurs accordingly condition is given by 


F C /jtanCo> - 1> + F^< /j + tanCa> > ■ 0 


DE: 

F 

« 0 

and F 

■= 0 ; 

X 

y 


EA: 

F 

= 0 

and u 

m 0 ; 


X y 


4St. STATIC ANALYSIS 

In this case initially, the die is completely flUed with 
the material and is assumed to be in the elastic state. Figure 4.4 
shows the biUet position in the die. Making use of the symmetry about 
the x-axis only half of the portion is analysed to reduce the 
computational effort. The mesh for this case is shown in figure 4.5 

The computational procedure is same as that of the 

transient case, except for the calculation of r^ 

In this case r Is calculated as follows 
3 

r m ^dlooxt (l^ . £) 

3 Au . 
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X 


Fig. 4.5 


Mesh for the F. E. M. 


analysis for static case 


3 




Where 


is the X" coordinate of the exit of the die 

dvaoxt 

x^ is the x-coordinate of the node which is nearer to the 
die exit and is on the die surface. 

All the above notations are shown in the figure ^iven bolow. 



Here r^ represents the factor that should give for the punch 
incremental displacement to Just move the i^^ node <whlch is on the 
surface of the die and nearer to the die exit> to the die exit. 

If r is minimum of r , r and r then i*^^ node shoxild move alon^ 

3 12 3 

DE as shown in figure 4.5 from the next iteration. 

BOUNDARY CONDITIONS 

AB: u ■ specified punch displacements and F * 0 ; 

X y 

BC: u * 0 and F ■■ 0 ; 

y 

CD: u tan<a> + u « 0 

X y 

and 

F C utanCa> -1> + F < /u+tan<ot> ) ■= 0 ; 

X ^ y 


DE; 

1 

o 

and 

F 

X 

« 0 ; 

EF; 

F « 0 

X 

and 

F 

y 

* 0 ■, 

FA: 

u » 0 
y 

and 

F 

X 

• *< 

0 

1 
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CHAPTER 5 


RESULTS AND DISCUSSIONS 

In -this chap'ter result.s for a V-no'tched problem are dlsctissed and 
compared wlbh -the Yamada-'sClI result/S -Lo vaUdabe -the pz*ocram and also 
results for the extrusion problem are presented. 

5.1 V-NOTGHED PROBLEM 

Fie- 5.1 gives the details of the v-notched tension specimen with 
a ■ 90° <notch an^^le). The material Is assumed to be elastlc-perfectly 
plastic and elastic linearly work hardening material. The material 
properties are as follows. 

Yoxm^^s modulus E • 20000 N/iran^ 

Poisson^s ratio i> ■ 0.3 

Yield sti^ess Y *30 N/mm* 

2 

S'train haixlenin^ paramnet^er H' * O N/'mm Cfor perfecrtly plasrtic> 

» 1000 K/mm^ bilinear maberial> 

P repr^sen't-si bhe load per half widt-h of bhe specimen. is bbe 

minimum cross secbional area of "the specimen. By making use of the 
symmetry about X*axls and Y-axis one quaKlrant of the sp>eclnien is 
anadysed to reduce the computational cost. The discretization of the 
specimen is shown in fig. 5.2 

Fig, 5.3 shows the development of the plastic zone at different 
stages of iterations. Numbers in the uppermost left corner represent 
the relevant stage number. 

Fig, 5.4 shows the variation of loadCPI> with respect to 
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deflect,lon<6> for the V-notched specimen. Deflection<<5> Is the 

deflection correspondinc the point C on the X-axis as shown in flp. 
S.2. It is observed that thet^ is no sudden chani^e in the overall 

stiffness of the specimen at the early stages. The load deflection 
curve has a sharp bend after certain sta^^e, when the plastic zone, 
spreading inward from both sides of the notch, meet first at a point 
on the longitudinal axis <X-axis>. After that plastic zone spreads 

very rapidly on both the sides. The specimen will fail at a stage, 
when all the element across the cross section Ci.e. the element along 
the Y-axis> aa>e yielded. 

Fig. 5.5 shows the distribution of the equivalent stress < a> 
contou 3 >s. From fig. 55 it is observed that the region near the root 
of the notch is highly stress concentrated. 

Fig. 5.6 shows the distribution of the longitudinal stressCo'^> 

over the minimum section. 

Fig. 5,7 shows the characteristics of the materials under 
consideration. 

The following section compaiies the results with YamadaCll. 
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Table 5.1 

Comparison of the results with those published by Yamada Cl] 



Present Method 

Y amada ' s Method 

Mac hi ne 

HP 9000 - 800 series 

HI TAG 5O20E 

No. of Nodes 

398 

' 144 

No. of elements 

735 

245 

Mesh pattern 

Hexagonal 

Hexagonal 

No. of stages 

43 

51 

Total computation 

22 min 

29 min 

ti me 



P 

234.6 N i 

221.643 N 




\ 

> 

0. 391 

0. 369 

e y 



P 

723. 316 N 

734.521 N 

c 



\ 

> 

1. 205 

1.224 

c y 



P / P 
c e 

3. 0897 

3. 314 
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5.2 EXTRUSION PROBLEMS 


The results for the following problems are presented in this 
secti on. 

<1> Transient state analysis of plane-strain extrusion problem 
<2> Transient state analysis of axisymmetri c extrusion problem 
<3> Static analysis of plane-strain extrusion problem 
<4> Static analysis of axi symmetric extrusion problem 
The prop>erties of the material under consideration are as follows 
y oung ’ s modul us E = 1 2500 N/mm 

Poisson's ratio ^ = 0.34 

Yield stress Y = 8 N/mm 

Strain hardening coefficient H' = O N/mm 

The deformed shapes for the different problems are shown in 
figures 5.8, 5.9 and 5.10. The dotted pattern shows the original shape 
of the material before going to deform. 

The plastic zone spreading "for the different problems at 
different stages are shown in figures 5.11, 5.12, 5.13 and 5.14. The 

numbers below the each of the figure indicates the number of elements 

in the plastic zone at that stage. 

The contour^ of the equivalent stress C & :> and equivalent strain 
CfiD for different problems are shown in figures 5.15 to 5.22. From the 
above figures it is observed that stress is highly concentrated at the 
die corners. 

And finally Fig. 5.23 shows the variations of total -work, 
elastic-work and plastic-work with displacement of the punch. From 
Fig. 5.23 it is observed that elastic work is not negligible compared 
to the total work as in the case of visco-plastic, model. 
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SPECIMEN 


romm ’ romm 



Fig. 5.1 V - notched tension specimen 
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Pig 5 2 Mesh for the V-notched problem 
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Fig. 5.5 Iso - eqxiivalent, stress contouiNS for V-notched 
specimen 
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Fig. 5.9 Deformed 
analysis 


mesh for transient. axisymmefric 
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Fig. 5.10 Deformed mesh for static plane-strain 
analysis 
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Fig. 5.11 Spreading of plastic zone tss- transient 
plane-strain case 
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I so— equivalent strain contours for transient 
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I^o-equi valent stress contours for transient 
a,>dL symmetric case 
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Fig. 5.18 I so -equivalent strain contours for transient 
axi symmetric case 
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Fig. 5.19 Iso-equivalent stress contours for static 
plane-strain case 
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WORK DISPLACEMENT 



Fig. 5.23 Graph between the differentjswork expenditure 
and the punch displacements 
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CHAPTER 6 


CONCLUSIONS AND SCOPE FOR FUTURE WORK 

CONCLUSIONS 

The following conclusions can be drawn from bhe present, sbudy. 

1. When t.he unloading Is not. considered at. t.he exit, of 'Lhe 
extrusion die, it is observed that the entire specimen yields. 

2. For conical die extrusion problem elastic work is comparad>le 
with plastic work. 

3. For axlsymmetrlc extrusion plastic zone is spreading more 
f aster than the plane strain extrusion at the exit of the die. 

SCOPE FOR FURTHER V/ORK 

1. Higher order elements may be used to take care of high stress 
gradient. 

2. Die geometry optimization. 

3. Because of plastic deformations, the variation in temperature 
and resulting change in properties may be taken care of. 

4. This Analysis may be extended to metal forming processes like 
forging, sheet metal forming etc. 
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5. Th« foUowini; Okreas of st.udy nuiy also be considered 


Hi^h stress concentration problems. 
Ductile fracture mechanics problems. 
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APPENDIX -I 


E)erivation for d\: 

Consider a simple uniaxial tensile test. 


d£* 


= dX 


Because 


Lc 


/ 


— A^. .Ar , 
3 


2 

3 


C A^^ - +C A£^ ~ At:^ C A^:^ - D' 

12 12 12 


Here Af^ = -2A£^^ = -2Ar^ 

12 3 


subscrip'ts 1 , 2 and 3 represen'Ls “the principal directions. 


As: ^ = Ae^ 


Substituting back to equation Cla!) 
We get , 


dc 


= dX 


d£: 


P _ 


dX cy' 

X 


= dX C cx - 


a + 0 * + O' 

X y 2 


Cla3 


] 


CIW 


CO (\] 


d\ 


a 

3 


X 


0 


d\ 


O' 


since a - o 
y z 


a 

3 


O' 


d.£ ^ ^ 3 do 

o So do/d£ 


dX 


1_ i. 

S - 


do 


H' 


CSD 


o 


H 



APPENDIX - II 


Equivalent stress a = 


^ J 


. + La\ .D Ccx' + La^ 5 

2 VJ LJ V.J \.j 


Where 


Acx is the increment of the stresses due to the increment of 


load. 


r Ao* = 

/I- 

Co: . + rAo: .5 

Co' . + rAo: .2) 


/ ^ 

tj i j 

J J 

is scaled 

by the 

factor r so 

that the equivalent 

is Just to 

reach yi 

eld stressCYD 


Y = + 

r Acx 



^ ^ \ c- 

o*' + a 

r O’' Ao*' . + 

r^ Ao: ^ ] 

2 L i j 

^ j 

i j 1- J 

vj J 

Aa^ 

— Aa' . 

Ao*' . 


st«p 

2 t J 

^ J 



— 2 — ■ ” 2 

o' +2r do AcyD 


ty ^ - Ac? f 

stop 


2 . “ 2 
+ r Ac? , 

stop 


Substituting for a + La from equation C2D , the above equation 
becomes 

-_2 “2 _ I- . ^ = A C3D 


Y^- cc^-arr + r^Ao”^ =0 

stop 
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Where 


r = Ao- * 2 <7 Aa - Lot * 

& t A p 


Equation C33 Is quadratic in r. By solving equation C3;5 
We get 


r 


r + 




8t«p 


Ae 


2 

stop 


C4D 
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